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[i]  This  paper  presents  a  convergence  study  for  second 
order  finite  difference  Laplacian  diffusion  used  in  ocean 
models.  For  demonstration,  ocean  model  simulations  are 
performed  over  a  rectangular  domain,  based  on  the  North 
Pacific  subtropical  gyre  region  with  grid  resolution  between 
1/2°  and  1/32°  and  with  horizontal  eddy  viscosity  coefficient 
(A//)  ranging  from  8000  to  30  m2  s  1 .  A  range  of  An  which  is 
appropriate  for  useful  model  simulations  of  an  oceanic 
domain  is  found  to  exist.  This  range  is  determined  by 
examining  the  spatial  patterns  of  Eddy  kinetic  energy  and 
mean  sea  surface  height.  The  results  fall  into  three  broad 
categories:  (a)  converged,  (b)  converging,  and  (c)  numerical 
problems.  Solutions  in  the  “converged”  category  do  not 
change  with  increased  grid  resolution,  and  solutions  in  the 
“numerical  problems”  category  exhibit  distinct  differences 
to  the  converged  result  at  the  same  A!h  Citation:  Wallcraft, 
A.  J.,  A.  B.  Kara,  and  H.  E.  Hurlburt  (2005),  Convergence  of 
Laplacian  diffusion  versus  resolution  of  an  ocean  model,  Geophys. 
Res.  Lett .,  32 ,  L07604,  doi:  1 0. 1 029/2005GL0225 1 4. 

1.  Introduction 

[2]  In  all  numerical  ocean  models  it  is  necessary  to 
parameterize  the  effects  of  unresolved  scales  of  motion. 
Since  the  horizontal  (isopycnal)  and  vertical  scales  are  so 
different  they  are  usually  treated  separately,  and  some 
form  of  Laplacian  diffusion  is  a  common  choice  for  the 
horizontal  parameterization  [e.g.,  Berloff  and  McWilliams , 
2002],  Diffusion  on  a  fixed  Eulerian  grid  is  complicated 
by  issues  concerning  the  difference  between  “horizontal” 
in  grid  coordinates  and  in  isopycnal  coordinates  [e.g., 
Kill  worth,  1997]. 

[1]  Layered  models  are  already  using  isopycnal,  or 
approximately  isopycnal,  coordinates  and  hence  give  a 
cleaner  separation  between  horizontal  and  vertical  diffusion 
[Smith  and  Gent ,  2004].  However,  Lagrangian  layers 
provide  several  possible  forms  for  the  Laplacian  diffusion 
term.  In  fact,  there  is  only  one  correct  form  which  is 
expressible  as  divergence  of  the  viscous  stress  tensor 
[e.g.,  Shchepetkin  and  O’Brien ,  1996].  In  Cartesian  coor¬ 
dinates,  the  form  of  Laplacian  diffusion  used  here  is  An  V2 
(h  V ),  and  the  correct  form  is  Afl  (V  ■  (hS7))  V ,  where  An  is 
the  horizontal  eddy  viscosity  coefficient  (m2  s“'),  V  is  the 
layer  velocity  (ms  ')  and  h  is  layer  thickness  (m).  These 
two  forms  are  equivalent  if  h  is  constant  as  would  be  the 
case  in  z-level  models  [e.g.,  Griffies  el  al 2000]. 

[4]  In  this  paper,  a  convergence  study  of  the  second  order 
finite  difference  horizontal  eddy  viscosity  variability  is 
undertaken  in  the  idealized  Pacific  Ocean  by  examining 
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simulations  from  an  ocean  model.  In  numerical  ocean 
modeling  studies,  “convergence”  can  have  at  least  two 
meanings.  Here,  we  have  taken  it  to  mean  convergence  to 
the  continuum  limit  of  the  finite  difference  equations  as 
the  horizontal  grid  resolution  is  refined.  An  alternative 
meaning,  of  more  practical  interest  to  ocean  modelers,  is 
convergence  on  the  actual  behavior  of  the  ocean  region 
under  study.  In  this  stronger  sense,  convergence  implies  no 
significant  difference  in  fidelity  to  the  real  ocean  between 
the  best  choice  of  model  parameters  at  one  resolution  and 
the  best  choice  of  model  parameters  at  a  finer  resolution. 
Here,  “best  choice”  can  obviously  include  a  lower  eddy 
viscosity  on  the  finer  grid  but  it  can  also  include  a  more 
accurate  coastline  and  bottom  topography. 

2.  Ocean  Model  Application 

[5]  The  numerical  model,  Naval  Research  Laboratory 
(NRL)  Layered  Ocean  Model  (NLOM),  uses  a  primitive 
equation  layered  formulation  where  the  equations  have 
been  vertically  integrated  through  each  Lagrangian  layer. 
Prognostic  variables  are  layer  density,  layer  thickness,  and 
layer  volume  transport  per  unit  width  (layer  velocity 
times  layer  thickness)  as  described  by  Wallcraft  et  al. 
[2003].  Our  experience  has  been  that  horizontal  diffusion 
is  largely  independent  of  the  number  of  layers  used.  Two 
isopycnal  layers  are  typically  sufficient  at  an  appropriate 
horizontal  resolution  and  value  of  An  to  yield  relevant 
(semi-quantitatively)  representations  of  mesoscale  eddies 
and  western  boundary  currents  in  subtropical-gyre-like 
model  domains  [ Schmitz  and  Thompson ,  1993].  The  sim¬ 
plicity  of  a  layered  model  formulation  permits  many 
numerical  experiments  to  be  performed  at  comparatively 
high  horizontal  resolutions  [ Hurlburt  and  Hogan ,  2000] 
with  a  variety  of  values  of  Aih  such  that  western  boundary 
current  path  and  structures,  including  recirculating  gyres  are 
also  realistically  represented.  Thus,  subtropical  gyre  simu¬ 
lations  in  this  study  use  a  two  layer  model,  a  configuration 
which  has  the  additional  advantage  that  mixing  between 
layers  (an  additional  source  of  damping)  is  minimized. 

[6]  Since  we  are  concerned  with  convergence  to  the 
continuum  limit  of  the  ocean  model  equations  in  this  study, 
the  model  domain  need  not  be  an  actual  ocean  basin. 
However,  the  results  must  be  representative  of  actual  ocean 
basins  since  these  are  our  primary  interest.  The  model 
simulations  are  performed  using  a  rectangular-shaped  area 
of  the  Pacific  Ocean,  corresponding  to  the  subtropical  gyre 
region,  with  closed  boundaries.  It  is  based  on  a  similar 
subtropical  gyre  region,  that  included  a  realistic  representa¬ 
tion  of  the  Japan/East  Sea.  For  this  study,  the  Japan/East  Sea 
has  been  converted  into  a  plateau.  Removing  the  Japan/East 
Sea  is  necessary  because  the  straits  are  not  resolvable  on  all 
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Table  1.  Grid  Resolutions  Used  in  the  Model  Simulations  of  a 
Rectangular  Pacific  Region3 


Resolution 

Ar  (deft) 

Ar  (deg) 

A.r  (km) 

Av  (km) 

1/2° 

0.703125 

0.5 

55.3 

55.6 

1/4° 

0.351569 

0.25 

27.6 

27.8 

1/8° 

0.175781 

0.125 

13.8 

13.9 

1/16° 

0.087890 

0.0625 

6.9 

6.9 

1/32° 

0.043945 

0.03125 

3.5 

3.5 

the  grids  used  here.  Rather  than  converting  Japan  to  a 
seamount,  it  would  have  been  possible  to  convert  the 
Japan/East  Sea  to  “land”.  However,  realistic  simulations 
of  the  Pacific  subtropical  gyre  are  not  possible  without  the 
Japan/East  Sea  [Hogan  and  Hurl  hurt,  2000].  Rectangular 
regions  are  less  expensive  to  run,  and  more  importantly  a 
rectangular  region  allows  prototyping  of  new  numerical 
schemes  without  requiring  implementation  of  general  coast¬ 
line  boundary  conditions.  The  disadvantage  of  a  rectangular 
region  is  that  it  may  increase  the  tendency  for  western 
boundary  currents  to  “overshoot”  the  expected  separation 
point  from  the  boundary.  However,  overshoot  is  also 
common  in  nonrectangular  regions. 

[7]  Here,  all  simulations  use  a  layered  Laplacian  diffu¬ 
sion  coefficient,  which  was  standard  in  earlier  versions  of 
the  model  [e.g.,  Hurlburt  et  al. ,  1996],  while  more  recent 
NLOM  simulations  of  actual  ocean  basins  use  stress  tensor 
Laplacian  diffusion  coefficient  [e.g.,  Kara  et  al .,  2003; 
Wallcraft  et  al .,  2003;  Kara  et  al .,  2004].  In  practice,  little 
difference  was  found  between  the  various  Laplacian  for¬ 
mulations.  In  particular,  repeating  a  few  simulations  from 
this  study  with  the  exact  form  of  Laplacian  diffusion 
showed  no  change  in  the  model  simulations.  A  layered 
Laplacian  diffusion  is  used  because  it  is  significantly  less 
computationally  expensive. 

[s]  The  model  simulations  are  performed  using  grid 
resolutions  of  1/2°,  1/4°,  1/8°,  1/16°  and  1/32°. 
Corresponding  grid  spacings  in  km  are  given  in  Table  1. 
Thirty  two  simulations  are  spun  up  to  statistical  equilibrium 
with  coefficients  of  horizontal  eddy  viscosity,  An ,  between 
8000  and  30  m2  s~ 1 .  Each  two  layer  simulation  uses 
identical  realistic  bottom  topography  from  a  modified  ver¬ 
sion  of  the  Earth  Topography  Five  Minute  Grid  (ETOP05) 
[National  Oceanic  and  Atmospheric  Administration ,  1988] 
and  the  wind  stress  forcing  is  Hellerman  and  Rosenstein 
[1983]. 

3.  Convergence  of  Eddy  Viscosity 

[9]  In  the  model  analysis  (see  section  4),  the  model 
simulations  are  subjectively  divided  into  three  broad 
categories:  (a)  converged,  (b)  converging,  and  (c)  numerical 
problems.  Since  simulations  are  nondeterministic  due  to 
flow  instabilities,  only  statistical  convergence  is  sought 
(i.e.,  mean  and  variability).  Solutions  in  the  “converged” 
category  do  not  change  statistically  with  increased  grid 
resolution.  They  are  essentially  identical  statistically  to  very 
high  resolution  at  the  same  eddy  viscosity.  Solutions  in  the 
“numerical  problems”  category  exhibit  distinct  differences 
to  the  converged  result  at  the  same  An.  Solutions  in  the 
“converging”  category  are  similar  to  the  converged  result, 
but  with  quantitative  differences.  It  is  possible  for  the 
converged  solution  to  change  drastically  as  the  eddy  vis¬ 
cosity  is  reduced.  This  can  happen  when  transitioning  from 


non-eddy-resolving  to  eddy-resolving  eddy  viscosities.  It 
can  also  happen  at  lower  eddy  viscosities  as  new  energy 
pathways  become  dominant. 

[10]  Figure  1  summarizes  the  results  on  a  graph  of 
Laplacian  diffusion  values  with  respect  to  grid  resolution. 
No  eddies  are  present  above  about  An  =  1000  m2  s  *,  and 
eddies  are  always  present  below  about  Au  =  500  m2  s  '. 
Since  the  categorization  is  subjective,  the  boundaries 
between  converged  and  converging  and  between  converg¬ 
ing  and  numerical  problems  are  not  sharply  defined  but  they 
appear  to  be  proportional  to  the  square  of  the  grid  resolution. 
The  open  circles  indicate  the  typical  eddy  viscosity  used, 
in  practice,  by  actual  NLOM  simulations  at  that  grid  reso¬ 
lution  [Hurlburt  and  Hogan ,  2000;  Hurlburt  et  al .,  1996]. 
At  1/2°  there  is  no  point  in  attempting  to  resolve  eddies,  and 
so  An  =  1 500  m2  s  1  is  typically  used.  At  1 14°  degree,  A/f- 
300  m2  s  1  is  very  marginally  eddy  resolving.  Marginally 
eddy  resolving  ocean  models  are  often  termed  eddy  permit¬ 
ting,  but  they  commonly  use  biharmonic  rather  than  La¬ 
placian  diffusion  [Gent  et  al. ,  2002;  Treguier  et  al .,  2001] 
and  are  therefore  outside  our  scope. 

[11]  Each  further  doubling  of  the  model  resolution 
decreases  the  typical  eddy  viscosity  by  a  factor  of  three 
(Figure  1).  All  the  typical  viscosities  are  in  the  “converg¬ 
ing”  zone,  because  we  are  attempting  to  “converge”  on  the 
actual  behavior  of  the  ocean  rather  than  the  continuum  limit 
of  the  finite  difference  model  equations.  As  grid  resolution 
increases,  the  typical  viscosity  value  used  in  the  NLOM  gets 
closer  to  the  converged  zone. 

[12]  In  practice,  for  regions  similar  to  that  used  in  this 
paper,  a  An  value  would  be  typically  chosen  in  the  converg¬ 
ing  category. 

4.  Eddy  Viscosity  Variations  With  Grid 
Resolution 

[13]  To  illustrate  the  results  given  in  Figure  1,  we 
examine  mean  sea  surface  height  (SSH),  a  good  indicator 
of  mean  flow,  and  eddy  kinetic  energy  (EK_E)  from  the  first 
model  layer,  a  good  indicator  of  flow  instabilities.  The 


Laplacian  diffusion  coefficient,  Alf  (m2  s  ') 

30  50  100  300  500  1  OCX)  2000  4000  8000 


Figure  1.  Laplacian  diffusion  coefficient  (also  known  as 
horizontal  eddy  viscosity)  values  with  respect  to  grid 
resolutions  used  for  the  OGCM  (NLOM)  simulations.  Solid 
dots  indicate  experiments  run  to  statistical  equilibrium,  and 
open  circles  indicate  the  “typical”  eddy  viscosity  that 
would  be  used  at  the  given  grid  resolution.  Note  that  the 
Laplacian  jliffiision  used  in  our  sequence  of  simulations  is 
An  V2(/?  V)  as  explained  in  the  text. 
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Figure  2.  Bottom  topography  (upper  left  panel)  for  the 
ocean  model  domain,  and  surface  layer  mean  Eddy  kinetic 
energy  (EKE)  per  unit  mass  fields  calculated  over  three 
years  from  the  model  simulations  configured  at  grid 
resolutions  of  1/2°,  1/4°  and  1/8°.  Mean  sea  surface  height 
(SSH)  is  overlain  on  EKE.  The  horizontal  eddy  viscosity 
coefficient  ( An )  value  of  2000  m2  s  1  is  used  in  the  ocean 
model  simulations.  Note  that  EKE  is  given  as  ‘dog  scale”  in 
the  color  bar.  SSH  is  plotted  using  solid  lines  (in  black)  with 
a  contour  interval  of  5  cm,  and  dotted  SSH  contours 
represent  negative  values,  i.e.,  less  than  the  domain  average. 
Long  dashed  contour  lines  are  used  for  zero  SSH.  Japan/ 
East  Sea  (JES)  has  been  converted  into  a  plateau  for  the 
model  simulations. 


means  of  SSH  and  EKE  are  calculated  over  three  years.  At 
least  a  three-year  mean  was  needed  because  the  finer 
resolution  NLOM  simulations  (e.g.,  finer  than  1/2°)  are 
nondeterministic  and  flow  instabilities  have  a  major  impact 
on  the  results.  The  simulated  SSH  and  EKE  fields  are  then 
compared  with  those  from  the  highest  grid  resolution 
available  at  the  same  eddy  viscosity. 

[u]  Figure  2  shows  EKE  and  SSH  at  three  different 
model  grid  resolutions  when  An  is  set  to  2000  m2  s  '.  The 
1/4°  and  1/8°  results  are  virtually  identical,  demonstrating 
convergence  at  1/4°.  The  1/2°  result  is  similar,  but  not 
identical,  to  those  at  higher  resolution;  indicating  that  this  is 
in  the  converging  category.  The  very  low  values  for  EKE 
are  typical  of  non-eddy-resolving  simulations.  When  An  — 
500  m2  s  1  is  used  in  the  simulations  at  four  grid  resolutions 
(Figure  3a),  the  results  from  the  1/8°  and  1/16°  model  are 
nearly  same,  indicating  an  almost  converged  solution  at 
1/8°.  Although  comparison  with  the  1/8°  simulation  indicates 
that  the  1/4°  simulation  has  serious  numerical  problems  near 
the  western  boundary,  we  place  it  in  the  converging  category 
(with  the  preceding  caveat)  because  1/4°  simulations  with 
realistic  geometry  for  this  region  do  not  exhibit  such  aberrant 
western  boundary  current  behavior,  just  qualitative  and 
quantitative  inaccuracy.  The  simulation  from  the  1/2°  model 
is  completely  different  again,  showing  the  narrow  zonal  jets 
that  are  often  associated  with  finite  difference  truncation  error 
and  is  therefore  placed  in  the  numerical  problems  category. 
However,  simulations  performed  with  An=  50  m2  s“ 1  clearly 
demonstrate  that  the  1/16°  and  1/32°  results  are  similar  but 
not  identical  (Figure  3b). 


[15]  No  1/64°  results  are  available  due  to  computational 
cost,  but  convergence  is  assumed  at  1/32°.  Despite  the 
severe  western  boundary  current  overshoot  in  the  1/8° 
simulation,  we  place  it  in  the  converging  category  for  the 
same  reason  as  the  1/4°  simulation  with  Alf  =  500  m2  s  1 
[see  Hurl  hurt  et  al. ,  1996,  Figure  7  and  Plate  8].  The  1/4° 
result  is  completely  different  again,  showing  almost  no 
penetration  of  the  primary  jet  into  the  basin  and  a  western 
boundary  overshoot  that  continues  along  the  northern 
boundary  and  part  of  the  eastern  boundary.  It  is  clearly  in 
the  numerical  problems  category. 

[16]  Comparing  the  last  panels  from  Figures  2,  3a,  and  3b 
demonstrates  the  large  impact  that  changing  the  eddy  viscos¬ 
ity  can  have  on  the  solution.  Suppose  the  1/32°  An  =  50  m2 
s" 1  solution  represents  the  real  state  of  an  actual  ocean  region. 
If  we  insist  on  running  simulations  that  are  fully  converged, 
then  at  1/8°  the  best  we  can  do  is  An  -  500  m2  s_l  which  is 
very  far  from  an  accurate  representation.  If  simulations 
are  allowed  in  the  converging  category,  then  the  1/8°  An  = 
50  m2  s-1  is  one  possibility  which  is  significantly  closer  than 
An  —  500  m2  s  1  to  the  desired  result.  In  practice,  we  would 
typically  use  An~  100  m2  s  1  at  1/8°  (not  shown),  because  it 
is  a  safer  distance  from  the  numerical  problem  zone. 

5.  Conclusions 

[17]  The  most  significant  results  from  this  paper  are: 
(I)  the  ocean  model  resolution  required  for  convergence 


Figure  3.  The  same  as  Figure  2  but  the  ocean  model  uses 
(a)  500  m2  s-1  and  (b)  50  m2  s' 1  at  grid  resolutions  ranging 
from  1/2°  to  1/32°. 
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of  Laplacian  diffusion  is  proportional  to  the  square  of  the 
grid  spacing;  (2)  it  is  not  safe  to  use  the  lowest  possible 
eddy  viscosity  because  there  is  high  risk  for  grossly 
inaccurate  solutions;  (3)  for  a  given  horizontal  model 
resolution  the  most  realistic  ocean  simulations  are  obtained 
when  the  eddy  viscosity  is  smaller  than  required  for 
convergence  to  the  continuum  limit  but  large  enough  to 
avoid  gross  finite  difference  truncation  error  effects  (the 
usual  practice  of  ocean  modelers),  and  (4)  the  typical  eddy 
viscosity  gets  closer  to  the  converged  value  as  resolution  is 
increased  (because  the  practical  eddy  viscosity  is  reduced 
by  a  factor  of  three,  not  four,  for  each  doubling  of  grid 
resolution). 

[ 1 8]  This  study  is  for  a  particular  model  in  a  particular 
idealized  region.  However,  our  experience  with  the  ocean 
model  used  in  this  paper  suggests  that  although  the  actual 
values  of  AH  used  can  be  region  dependent  (based  on 
maximum  current  speeds),  the  results  are  insensitive  to 
vertical  resolution,  and  that  the  same  general  patterns 
illustrated  here  apply  to  all  ocean  regions.  NLOM  has  fewer 
sources  of  diffusion  than  some  other  ocean  models,  but 
these  results  are  robust  enough  to  suggest  that  they 
may  apply  to  most  if  not  all  ocean  models  that  combine 
Laplacian  diffusion  with  second  order  finite  differences  in 
space. 
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